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Abstract 

We analyze a biophysical model of a neuron from the entorhinal cortex that includes persistent 

sodium and slow potassium as non-standard currents using reduction of dimension and dynamical 

, systems techniques to determine the mechanisms for the generation of mixed-mode oscillations. We 

. have found that the standard spiking currents (sodium and potassium) play a critical role in the 

OO ■ analysis of the interspike interval. To study the mixed-mode oscillations, the six dimensional model 

O ' 

• ' has been reduced to a three dimensional model for the subthreshold regime. Additional transforma- 

' tions and a truncation have led to a simplified model system with three timescales that retains many 

00 ! 

, properties of the original equations, and we employ this system to elucidate the underlying structure 

and explain a novel mechanism for the generation of mixed-mode oscillations based on the canard 
^ . phenomenon. In particular, we prove the existence of a special solution, a singular primary canard, 

■ that serves as a transition between mixed-mode oscillations and spiking in the singular limit by em- 

ploying appropriate rescalings, center manifold reductions, and energy arguments. Additionally, we 
conjecture that the singular canard solution is the limit of a family of canards and provide numerical 
evidence for the conjecture. 
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In this paper we investigate a specific biophysical model of a neuron, not pre- 
viously analyzed, using techniques from dynamical systems theory. We discover 
the presence of mixed-mode oscillations (MMOs), alternating sequences of large 
amplitude (neuronal spikes) and small amplitude (subthreshold) oscillations, in 
the model. To explain the dynamical mechanisms responsible for the MMOs, we 
apply reduction of dimensions techniques to reduce the original model to a sim- 
plified model system with three timescales. Our analysis of this system leads to 
a novel mechanism for the generation of MMOs that is based on the presence of 
canard solutions that act as a transition between MMOs and spiking. The theory 
that we develop may be applied to systems from various contexts that exhibit 
mixed-mode oscillations and have the same underlying dynamical structure. 



INTRODUCTION 



Mixed-mode oscillations (MMOs) refer to oscillatory behavior in which a number (L) of 
large amplitude oscillations is followed by a number (s) of small amplitude oscillations, differing 
in an order of magnitude, to form a complex pattern (Lj^Lg^ . . . Lf^) of oscillations of varying 
amplitudes. MMOs were first noticed in the Belousov-Zhabotinsky reaction [l| , and have 
been subsequently witnessed in experiments and models of numerous chemical HQ, H, E| and 
biological Ja ^, a| systems. Recently there has been work on MMOs in models of neurons 
121 . Il3l . Il4| . In the context of neural dynamics, large oscillations correspond to firing 




11 



of action potentials while small oscillations correspond to subthreshold oscillations (STO). 

We have discovered MMOs in a biophysical model of a neuron in the entorhinal cortex 
(EC) introduced by Acker et al. 15||. The EC is thought to play a prominent role in memory 
formation through its reciprocal interactions with the hippocampus and participation in the 
generation of the theta rhythm (4-12 Hz) 16, l3, 18, 19]. Pyramidal neurons in layer V of the 



entorhinal cortex have been found to possess electrophysiological characteristics very similar to 



those seen in layer II stellate cells of the entorhinal cortex |20|, |21[, including their propensity 
to produce intrinsic subthreshold oscillations at theta frequencies. In layer V pyramidal cells 
these subthreshold oscillations are thought to arise from an interplay between a persistent 
sodium current and a slow potassium current 0, 0], while in layer II stellate cells, they are 
thought to arise from an interplay between a persistent sodium current and a hyperpolarization 
activated inward current (Ih) [1 



3, El- 



The model from Acker et al. ISj that we employ includes these currents and captures some 
dynamical aspects of the generation of STOs in layer V pyramidal cells. After a numerical 
investigation of the dynamics of the full model that displays the interaction of the persistent 
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sodium and slow potassium currents to form STOs, we modify it by replacing the persistent 
sodium current by a tonic applied current. This can be interpreted as the modeling of an 
electrophysiological experiment designed to investigate the effects of blocking the persistent 
sodium current. The modified model displays MMOs, and the main focus of this work is to 
understand their nature. In order to simplify the analysis of the model, we perform a reduction 



of dimensions analysis 12|, |2J] to reduce the 6-d system into a 3-d system. 



Among the proposed mechanisms for MMOs are break-up of an invariant torus [2q] , break 



up/loss of stability of a Shilnikov homoclinic or bitp, l26l | and subcritical Hopf bifurcation 



combined with a suitable return mechanism 27|, |28|]. Well known in the context systems 



itji, 
Si. 



with multiple time scales is the canard mechanism, introduced (to the best of our knowledge) 
by Milik and Szmolyan j^. It is based on the presence of canards, which are defined as 
solutions that originate in a stable slow manifold and continue to an unstable slow manifold 
while passing through a thin region of non-hyperbolic behavior, which we refer to as the fold 
region. Primary canard solutions, that neither perform small amplitude oscillations nor spike, 
can act as dynamic boundaries between regions of small oscillations and spiking. Canard 
solutions have been analyzed in three dimensions using non-standard analysis by Benoit 29 1 
and using techniques of geometric singular perturbation theory in combination with blow-up 



techniques by Szmolyan and Wechselberger 30|, l31[. A general mechanism of MMOs based on 



the existence of canards has been postulated in 32] ■ There is a variety of different examples 



which, with some modifications, fit into the context of this general mechanism 

3, 0,13,0, Q- 

Our model fits into the context of a system with multiple timescales; we will show that 
the occurrence of MMOs in our model is related to the presence of canards. However, while 
in most of the other examples of canard mechanism small oscillations are due to rotations 
about a canard solution called the weak canard 32], in our case, during the stage of small 
oscillations, some of the trajectories follow closely the unstable manifold of an equilibrium of 
oscillatory type, very close to a Hopf bifurcation. In this aspect our system is similar to the 



systems studied in j2,|27|, 128|. Consequently our model provides a novel mechanism for MMOs, 
one in which both canards and dynamics near a saddle-focus equilibrium play an important 
role. 

We designed a three time scale model system that accurately reproduces most of the fea- 
tures of the neuronal model, and a significant part of this work is devoted to the analysis of 
the three time scale model equation. Most of the work, at least to some degree, relies on nu- 
merics. We show that canards can exist, but we conjecture they occur on parameter intervals 
of exponentially small width {0{e~^^'^), where e is the ratio of the two slowest timescales). 
This is similar to the case of two dimensional canard explosion (known also as the canard 
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phenomenon), where a small oscillation grows through a sequence of canard cycles to a re- 
laxation oscillation as the control parameter moves through an interval of exponentially small 
width 34, 3^, 36|. In contrast, in the two slow/one fast variable systems studied by Szmolyan, 



Wechselberger and others, canards occur in a robust manner, i.e., on open parameter intervals 
of size independent of the singular parameter. We also discover the role of canards in our 
system: we show (relying to some degree on numerics) that the transition from MMOs to 
pure spiking must correspond to a canard solution. We compute a canard numerically and 
compare the relevant parameter value to the parameter value of the transition from MMOs to 
spiking, finding that the two values are almost equal. We also calculate analytically the e = 
approximation of the numerically computed canard. 

An outline of the paper is the following. In section [Tll we introduce the model and display 
the dynamical properties and bifurcation structure of the original Acker et al. model through 
simulations of the full system. We also perform a reduction of dimensions analysis to reduce 
the 6-d system to a 3-d system with three time scales in the subthreshold regime. In section UTTt 
we introduce the modification of replacing the persistent sodium current with a tonic applied 
current, and we analyze the mechanism for the MMOs observed in the modified model. We 
begin with a phase space description of the MMO activity through simulations. Next we 
propose a three time scale model problem, by adding an additional slower time scale to the 
model of Szmolyan and Wechselberger 30|], and transform our system into this form. The 
three time scale model has a Hopf bifurcation and some of the dynamics in MMO sequences 
follows closely the unstable manifold of the unstable equilibrium, which reproduces well the 
behavior of the neuronal model. In our analysis of the MMOs in the three time scale model, 
a fundamental observation is that during the approach to the fold region the trajectories are 
attracted to a single trajectory, a 1-d (super) slow manifold contained in the attracting branch 
of a folded 2-d slow manifold. To investigate the flow near the fold curve, we perform an e 
(singular perturbation parameter) dependent rescaling which blows up the flow near the fold 
curve. We define singular canards as solutions of this system satisfying certain asymptotic 
properties and find explicitly a singular canard that is a polynomial vector function of the 
independent variable. By employing appropriate rescalings, center manifold reductions, and 
energy arguments, we prove that a continuous family of canard solutions must converge, as 
e ^ 0, to a singular canard. Consequently, we prove that if a transition from MMO to spiking 
occurs along a continuous curve ^{e) then /i(0) must be a value for which a singular canard 
exists. We conjecture that the mentioned polynomial solution is the only existing singular 
canard and provide numerical evidence for this conjecture. We conclude with a discussion of 
related and future work in section [IVl 
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II. THE MODEL 



A. Model equations 



To investigate STOs and MMOs in EC layer V pyramidal cells, we consider a biophysical 
model introduced by Acker et al. [15| that includes two currents (persistent sodium (iNap) 



and slow, non-inactivating potassium (/ 



the generation of STOs in t 
Hodgkin-Huxley currents 



lese cells 
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whose interaction is known to be responsible for 
8|l- In addition, the model contains the standard 



: transient sodium (/atq), delayed rectifier potassium (Ir), and 



leak (/l). The current balance equation for the model is given by 



C— = -QNam^Kv - ENa) - gKn^{v - Ek) - giiv - El) (1) 

-9NapP{v - Eno) - 9KsW{v - Ek) + hpp 

where C is the membrane capacitance {fiF/cm^), v is the membrane potential (mV), lapp is 
the applied bias current (fiA/cm?), g is conductance {mS/crn?'), E is the reversal potential 
(mV), and the time t is in ms. Also, we define /TVa = QNa'i^'^hiy — Ej^a), Ik = QK^'^iy — Ek), 
h = giiv - El), iNap = gNapP{v - Emcl), Irs = OksW^v - Er)- The values of the parameters 
can be found in the Appendix. The dynamics of each of the gating variables m, h, n,p, and w 
is described by a first order equation of the form 

dx Xqq(^v) X 

dt T^{v) 

where x represents the gating variable, Tx{v) represents the volt age- dependent time constant, 
and Xoo{v) represents the voltage dependent steady state values of the gating variable. We will 
see that the iNap and Iks play an important role in the generation of subthreshold oscillations. 
The activation and inactivation curves Xoo{v) for the gating variables are shown in Figured] 
while the voltage dependent timescales Tx{v) are shown in Figure [2l 



B. Simulations of the full system 

Simulations of the full system ([I])-® show that the model fits into the category of class 2 
neural excitability js^ . Using the applied current as a bifurcation parameter, the system un- 
dergoes a subcritical Hopf bifurcation as lapp is increased as seen by Figure [31 Thus, as lapp is 
increased, the cell switches between three different regimes . In the first regime {lapp < 0.88), 
there is a single stable rest state which is converged to through damped subthreshold oscilla- 
tions. In the second regime (0.88 < lapp < 0.994), there is bi-stability of the rest state and 
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periodically spiking solutions, which are separated by unstable, small amplitude oscillations 
generated through a subcritical Hopf bifurcation (See Figure H]). The spiking solutions start 
with a minimal frequency of about 5Hz while the damped subthreshold oscillations have a 
frequency of about lOHz. Also, the replacement of the slow potassium current with a tonic 
drive results in the absence of subthreshold oscillations (as in j^) while the replacement of 
the persistent sodium current with a tonic drive will be discussed in section IIIII We note 
that the damped subthreshold oscillations involve the interplay of the slow potassium current 
and the persistent sodium current as seen in Figure [5l At the transition between the first 
and second regimes, we observe dynamics reminiscent of a 2-d canard phenomenon in which 
the trajectory shown follows a repelling unstable manifold. The third regime {lapp > 0.994) 
consists of stable, periodic firing. 

C. Reduction of Dimensions 

We show that during the interspike interval, which corresponds to the subthreshold voltage 
regime, we can reduce the system to a 3-dimensional system ([3]). By analyzing the voltage 
dependent timescales r^i^v) seen in Figure[2], we see that r^i^v) » T„(f ), Th{v) >> rm{v), Tp{y). 
Thus, we seem to have three timescales with p and m evolving very quickly and w evolving 
very slowly. Thus, letting p = Poo{v) and m = m^{v) yields a reasonable approximation. By 
considering Figure [1] and [2], we see that when v is hyperpolarized, h evolves much faster to its 
steady state value (its time scale is 4 times smaller) than n. Also, we note that n'^{v) > h'^{v). 
Thus, h quickly approaches its steady state value hoo{v), and then it can track it closely since 
hoo{v) changes very slowly during the interspike interval. Thus, we replace h with hoo{v) and 
retain n. Hence, we can reduce the system to one for which only v, w, and n are dynamic 
variables. 

To investigate the importance of the spiking currents (sodium and potassium), we also 
perform a bifurcation analysis. When /jv^ is removed, the system exhibits a supercritical Hopf 
bifurcation in which sustained, stable subthreshold oscillations appear. (See Figure O) This 
change in the bifurcation criticality due to Ing removal shows that the transient sodium current 
plays an important role in the interspike interval. With Ik removed, the bifurcation structure 
is retained, but there is a significant shift in the Hopf bifurcation point (from lapp ~ .9941 to 
lapp ~ .775). By retaining the n equation, and replacing m, h, and p with their steady state 
activation and inactivation functions moo{v) , hoo{v) , and Poo{v), respectively, we obtain the 
reduced 3d system 
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dt 



dw 
'dt 
dn 
Itt 



-9Naml^{v)hoo{v){v - Eno) - 9Kn^{v ■ 

-9NapPoo{v){v - En a) - gKsW{v - Ek) 



- Ek 

+ 1 



9l[v 



Er 



app 



(3) 



W 



n 



where Tw{v) has been replaced by its constant value, r^. We note that the bifurcation diagram 
for this system (see Figure [7j) is very similar to that of the full system. 

Typically, the spiking currents are not included in the biophysical description of the sub- 
threshold regime as in generalized integrate-and-fire models, which are constructed by bio- 
physically describing the subthreshold regime using non-spiking currents and adding artificial 



spikes 12|, |39|]. This system provides an example where neglecting the spiking currents in the 



biophysical description of the subthreshold regime leads to a qualitative change in the cell's 
dynamics, i.e. Ik and I^a must be included in this reduced system since our model is sensitive 
to small changes in the sodium and potassium currents during the latter part of the interspike 
interval before the spike. In some cases it has been shown that these currents play no critical 
role in the subthreshold regime 121]. 



III. ANALYSIS OF MIXED-MODE OSCILLATIONS 
A. Simulations and outline of the results 

We have discovered that in the full system and in the reduced system the re- 

placement of the persistent sodium current with a tonic excitatory drive leads to mixed-mode 
oscillations in which a large amplitude spike is followed by a number of small amplitude sub- 
threshold oscillations (see Figure El). We note that this modified version of ([3]), with the 
persistent sodium current removed, exhibits a more complex bifurcation structure than ([3]) it- 
self as seen in Figure [91 The system exhibits a single equilibrium which loses stability through 
a subcritical Hopf bifurcation at lapp ~ 16.93. For lapp ~ 18.25, the system exhibits continuous 
spiking. For lapp between ^ 17.61 and ~ 18.25, the system displays MMOs of type 1^. For 
lapp between ^ 16.98 and ~ 17.61, the system displays MMOs with multiple STOs per spike, 
and the number of STOs per spike is a decreasing function of lapp- Finally, for lapp between 
^ 16.9 and ~ 16.98, the system exhibits sustained subthreshold oscillations. 

Our goal is to analyze the mechanisms for this behavior in with qna = 0. We begin 
with a geometric interpretation of the MMOs using simulations of (|3]). We note that the 
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v-nullsurface is a 2-d surface w{v, n), and for fixed n, w{v, n) is a cubic sliaped curve. Also, we 
note that ([3]) possesses a single fixed point, denoted by Pq = {vo,Wo,nQ), which is near the fold 
curve of the v-nullsurface. The linearization about Pq yields a large negative eigenvalue and a 
pair of complex eigenvalues with positive real part, and the ratio of the real eigenvalue to the 
real part of the complex eigenvalue is about -30. Thus, Pq possesses a 1-d stable manifold and 
a 2-d unstable manifold. Numerical simulations show that a MMO trajectory travels down 
near the left sheet of the v-nullsurface (slow manifold) and subsequently passes near the fold, 
where it gets attracted to the 2-d unstable manifold of the fixed point Pq. See Figure fTOl The 
1-d stable eigenspace of Pq is plotted in red in the figures, and the plane displayed in the figures 
is the unstable eigenspace of Pq. The trajectory approaches the unstable manifold of Pq, which 
we denote by Mq, along the stable manifold of Pq. Then the trajectory spirals out from near Pq 
along Mq and continues to transversally pierce the v-nullsurface (the loss of alignment with the 
fast nuUsurface during the passage through the fold region is typical for folded node or folded 
saddle node type dynamics), until it becomes tangent to the v-nullsurface. At this point, the 
trajectory leaves Mq and gets attracted to the right sheet of the v-nullsurface since it enters a 
region in which v is strongly repelling from the left sheet. Then, there is a return mechanism 
(spiking regime) that brings the trajectories back to the left sheet of the v-nullsurface. We 
note that the feature of a passage near the unstable manifold of a saddle focus equilibrium 
followed by a reinjection mechanism is reminiscent of Shilnikov's example 401 as well as of 
the examples where MMOs occur through a subcritical Hopf bifurcation (see section 

|lV]for a detailed discussion). 

In the remaining part of this work much of the focus will be on the following model system 

e'^x = —y + + Dix^ 

y = X — z (4) 
ez = fi + ax + bz, 

where e > is small and fi is a. control parameter. System (jl]) shares many of the features 
of ([3]) but is a lot simpler. The initial part of our analysis is to establish two basic geometric 
features of the flow: 

1. There exists an equilibrium of oscillatory type, which loses stability through a Hopf 
bifurcation for close to 0. 

2. The critical manifold y = x"^ + Dix"^ is a surface with two fold lines. There exists a 2D 
attracting slow manifold close to one of the sheets of this surface, and within this slow 
manifold there is a ID super-slow manifold, approximating the curve y = x'^ + Dix^, 
X = z. 
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We will show that all the MMO trajectories must pass very close to the super-slow manifold. 
Consequently, the distance between the continuation of the super-slow manifold to the fold 
region and the stable manifold of the equilibrium is of great importance for the nature of the 
MMOs. If this distance is small, then the trajectory must follow for a long time the flow along 
the 2D unstable manifold of the equilibrium, making many rotations about the equilibrium. 
As this distance increases, the role of the unstable manifold of the equilibrium diminishes. 

In this work we have focused on the transition from MMOs to spiking. This occurs in 
the parameter regime where the continuation of the super-slow manifold, and thus the MMO 
trajectories, are relatively far from the equilibrium. We have been able to show, with the help 
of some numerics, that the transition from MMOs to spiking occurs through a canard, and 
that families of canard solutions converge to singular canards, which are special solutions of 
a suitable limiting system. We have also found a parameter value for which a special solution 
exists. We conjecture, based on numerical evidence, that this parameter value is unique. 



B. Derivation of the model system 

In order to investigate MMOs in system ([3]), we transform it to a three timescale toy model 
of the form 

x' = —y + F{x) 

y' = e\x - z) (5) 
z = e{fi + ax + hz) 

with F{x) = x"^ + Dix^ and e ^ 1. To this end, we begin by first treating ([3]) as a slow-fast 
system with two timescales in order to obtain a folded singularity, which we will use as the 
origin for our model system in order to study small amplitude oscillations near the fold curve 
of the 2-d folded v-nuUsurface of ([3]). For an exposition on folded singularities, see [30|, ISlI ]. 

Two timescales are uncovered by dividing the v equation of ([3]) by the maximum conduc- 
tance Qiqa and setting k = C/ g^a obtain 

K— = -ml^{v)h^{v){v - E^a) -QKn^iv - Ek) -Qhiv - El) 

-9NapPoo{v){v - Eno) - 9KsW{v - Ek) + lapp = /{v, W, u) (6) 

dw Woo{v) — w 



dt Tw 
dn n^iy) — n 

dt Tn{v) 



= 92iv,w,n) 
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where = Qx/QNa for x = K,L,Nap, and Ks and /, 



app 



Na- 



We note that the 



conductances in the v equation are bounded by 1 and the magnitudes of the w and n equations 
are bounded by 1 since = 90 and 1 < r„(i;) < 6. Thus, v evolves on a faster time scale 
then w or n. Momentarily, we neglect the fact that w and n vary on different timescales. 



By setting k = 0, we obtain the reduced system which is a 2-d dynamical 
critical manifold S = {{v,w,n) : f{v,w,n) = 0}. Based on the work of 30|, 



system on the 



31| . we wish to 



show that 5 is a non-degenerate folded surface in the neighborhood of a folded equilibrium 



point {vc,Wc,nc) with f{vc 



0. 



0, fy{Vc, Wc, Tie) = 0, fw{Vc, W, 

In order to compute {VcWcUc), we begin by solving f{v,w,n) = for w and obtain w = 
W{v,n). Next, we solve fv{v,W{v,n),n) = for f as a function of n and find that f is a 
constant, which we denote by Vc, along the fold curve (since / and are linear in n'^ and w). 
Also, we denote W{vc, n) by ^{n). Thus, the fold curve is parametrized by (t>c, ^(n), n), n E I. 
Also, we differentiate f{v,w,n) = with respect to t to obtain f^v = —fwQi — /n5'2- Thus, 
the reduced 2-d system on the critical manifold can be denoted by 



fvV = fwQl + fng2 

n = 92 



(7) 



where w is replaced by W{v,n). Through a solution dependent time rescaling we obtain the 
desingularized reduced system 



V = fu,9l + fn92 

n = -fv92- 



(8) 



We numerically solve fw9i + fn92 = for n with v = Vc and w = $,{n) to obtain Hc- Finally, 
"W^c = ^{nc). We note that {vc,Wc,nc) is a fixed point of ([7]) and ([8]). We linearize ([8]) about 
this fixed point and find that in our parameter regime both eigenvalues are negative, resulting 
in a folded node equilibrium. We note that the quotient of these eigenvalues predicts the max- 
imum number of subthreshold oscillations per spike 30, El| for k sufficiently small. However, 
this prediction is incorrect (underestimate) for this system since there are three timescales, 
resulting in a k that is not sufficiently small in the two timescale formulation. 

We now shift the folded node point to the origin in ([6]) with the transformation v = v — Vc, 
w = w — Wc, and n = n — ric to obtain 



Kv = f{v + Vc,w + Wc,n + nc) = f{v,w,n) 
w = 9i{v + Vc,w + Wc,n + nc) = gi{v,w,n) 
n = g2{v + Vc,w + Wc, n + Uc) = g2{v, w, n). 



(9) 
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Next, we rectify the fold curve to the n axis. We note that the fold curve is parametrized 
by (0, where x{^) = W{vc,fi + ric) — Wc- Thus, we make the substitution v = 

w = w — and ri = n to obtain 

Kv = f{v,w + xin),n) = f{v,w,n) (10) 

We note that /(O) = = #i(0). 

Then we let v = jiv, w = 72W, and n = n where 71 = |/^t,(0) and 72 = ^/to(0)/^(0) to 
obtain 

Kv = V + w'^ + 0{vw,v'^) = f{v,w,n) (11) 
n = g2{v,w,n). 

To simplify the n equation we let f = -u, w = w, and n = n/g2{Q) to obtain 

Kv = v'^ + w + 0{vw,v'^) 

w = [3n — cv + Oiw^n^ ^nv.v'^) (12) 
h = 1 + dn + ev + 0{w,v'^,nv). 

Next, we rescale time with t = t/ n, let x = v,y = —w, and z = {f3/c)n, and truncate higher 
order terms to obtain 

x' = -y + x^ + Dix^ 

y' = Kc{x — z) (13) 
z = —(1 + + ex). 

Finally, if we let e = ^/kc, then we obtain 

x' = -y + x'^ + Dix^ 

y' = e^{x-z) (14) 
z' = e(/i + ax + hz) 

where /i = /^^kc"^/^), a = /Jey^c"'^/^, and b = d^/nc^^^'^. This system is equivalent to ([5]). 
For example, if lapp = 17.1 in ([3]), the parameters are approximated by: e = .2764436178 = 
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Ed, fi = .02839202004, a = .4266759683, and b = -.9420074624. Also, we consider = -A, 
which provides an appropriate return mechanism. A more general assumption that we make 
to guarantee the desired functioning of the toy model is 

-fe>a>0. (15) 



C. Slow manifolds of the model system 



Equation ([5]) has a two dimensional critical manifold defined as the set of equilibria of ([5]) 
with e set to 0. We denote this manifold by S and note that it is defined by the equation 
y = x'^ + Dix^. Let x^ax denote the unique value of x for which the function F{x) = x'^ + Dix^ 
has a maximum. Let 



Sa = {{x,y, z) e S : X <0} 
Sr = {{x,y, z) E S : < x < x„ 



where a and r denote attracting and repelling, respectively. By Fenichel theory j41(], away 



from the fold lines x 



and X 



y = fi^max), system ([5]) with e > has locally 
invariant manifolds Sa,£ and Sr^e that are 0{e) close to Sa and Sr- While the existence of 
Sa^s and Sr^s is standard and not surprising, we also argue that there exists a one-dimensional 
super-slow manifold (it can be chosen to be contained in Sa^e), denoted by SiD,e, that attracts 
all the dynamics on Sa,e- (See Figure [TTl) For this reason it is convenient to rescale ([5]) to the 
slow time scale, which gives (jl]). Setting e = 0, we obtain two constraints: 



-y + x'^ + Dix^ 
fi + ax + bz 

This gives a one dimensional curve defined by: 








y 



x^ + Dxx^ 
fi + ax 



(16) 



which we denote by Sid- Note that assumption ( fT5l) guarantees that the flow on Sm is towards 
the fold. 

Linearizing (jl]) about Sid we obtain the matrix 

/ i(2x + 3Dia;2) 
A= 1 



\ 



-4 \ 
-1 
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Ignoring the 0(1) terms in A we obtain 

/ ^(2x + 3Dix2) 0\ 
0. 

V ! of/ 

It is easy to verify that this matrix, apart from an eigenvalue 0, has eigenvalues ^(2x-|-3DiX^) 
and |, which are both negative for a; < 0. It follows that the dominant eigenvalues of A are 

ir(2x + 3Dix'^) + 0(l) and - + 0(1). 

e 

There is also a third eigenvalue, which is 0(1), but we are interested only in the two leading 
ones. By a slight modification of Fenichel theory we can prove that ^ has a one dimensional 
attracting locally invariant manifold, which we denote by SiD,e- All the trajectories which 
make a relaxation loop (spike) must become exponentially close to it before they enter the 
fold region. Thus, there is a single trajectory that needs to be tracked through the fold region 
to determine the existence of MMOs/small amplitude dynamics and this trajectory can be 
arbitrarily chosen as long as its initial condition is not too close to the fold. We illustrate this 
scenario in Figure [12] in which we display Sir, along with a spiking solution for e = .1 that 
tracks Sid quite closely as well as a MMO solution for the larger value oi e = Sd- 



D. The microscope: a rescaling of the model system 



In order to understand the flow of ([5]) in the fold region, we introduce a rescaling which 
blows up a rectangle of size 0{e) x 0(6^) x 0{e) x 0(e) to the size of 0(1) in each of the 
directions. This trick has been widely used in the study of canard problems and was given the 
name 'microscope' by Benoit, Diener and others who investigated canard explosion and related 
problems by means of non-standard analysis 29|, 
rescaling is as follows: 

X = ex, y = e'^y, z = ez, jj = eft. 



431]. In the case of (JSj), the 'microscope' 



Remark 1 In Proposition [T] we will prove that MMOs cannot exist if is too large. This 
implies that MMOs exist for fi = 0{e) only. 

After rescaling the variables we obtain the following system: 

x' = - y + x'^ + eDix^ 

y' =x — z (17) 
z' =fi + ax + bz. 
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By setting e = we obtain: 



X 



— y + X 



y' =x — z 

z' =fi + ax + bz. 



Remark 2 A similar normal form arises in the study of folded node 29j, 30|]. The difference 
is in the z' equation, where the terms ax and bz do not appear. However these two terms 
alter the dynamics dramatically, so that the results on the folded node normal form cannot be 
easily applied in the context of flTSl) . We have some partial results in the direction of reducing 
the study of the dynamics of flTSl) to facts about the dynamics of the normal form of folded 
node, but the arguments are very technical and it does not seem worthwhile to present these 
results in this paper. 



The system ( !T8|) has a unique equilibrium given by 



X 



y 



a + b \a + b 

The linearization at the equilibrium is given as follows: 

_1 \ 



An 



\ 



2fi 

" a+b 

1 

a 








-1 



with det Aq = a + b. We assume that (a + b) < 0, which implies that the determinant is 
always negative. Clearly the only way for the matrix to be non-hyperbolic is if there is a pure 
imaginary eigenvalue. Let p{X) denote the characteristic polynomial of Aq and note that p{X) 
can be written in the form p{X) = A'^ + 02-^^ + ^lA + a^. The condition for Aq to have a purely 
imaginary eigenvalue is ao = 0102. We evaluate the coefficients and get 



ao 



(a + b), ai 



-2fLb + a + b 



a + b 



02 



-b{a + 6) + 2/2 
a + b 



We define the function C(/i) by the formula 



C{fi) = (ao — aia2){a + b). 
Clearly the roots of C(/i) give purely imaginary eigenvalues. Evaluating C(/i) we get 

C{ji) = {a + b)a + (26^ + 2)Ji 



46 



a + b 
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Now C(/i) has two roots, given by 

(1 + 6^ ± ^(l + b^y + Aba)^. 

For example, if a = .4266759683 and b = -.9420074624 (corresponding to hpp = 17.1), the 
two roots are fl = 0.06692624850 and /2 = 0.4493251582. The first corresponds to the Hopf 
bifurcation point while the second is an extraneous solution. 



E. Special solution 



System (jlSil has special solutions that are algebraic as t — *^ —oo, i.e. solutions that tend 
to — oo no faster than a polynomial function of t. In this section we will prove that for each 
value of p, there exists exactly one such solution (up to time rescaling). The importance of 
special solutions will become clear in Section IIII Fl where we will show that the continuation 
of the manifold SiD,e to the region of validity of the 'microscope coordinates' {x,y,z) is 0{e) 
close to the special solution of ( |T8l) corresponding to the same value of fl. 

In order to understand the flow of (|T8|l near infinity, including special solutions, we make 
the following transformation: 

cr = —T=, X = ax, z = az. (19) 
The transformation yields the following system: 

X = a (~2^ ~ '^'^ — 1 + X ) 
1 

a' = — a~(x — z 
2 ^ 

z! = aji \ ax ■\- bz — -^(^{i — z)z. 

Through an appropriate (solution dependent) time rescaling we can eliminate a factor of 
in ([2U]): 

x' = — cr'^ix — z)x — 1 + 
2 ^ ' 

1 



'- "^2,^ (20) 



cr 



2 



a\i-i) (21) 



/ = (T((T/i + ax + bz ~ -cr(x — z)z). 
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Solutions of fl?T]) . suitably reparametrized, become solutions of fl2U]) . except of course for the 
solutions for which cr = 0. Equation (|2T1) has two lines of equilibria, given by 



x = ±l, a = 0. (22) 

These lines of equilibria were artificially introduced by the time rescaling which led from fl20l) 
to (12T]) . The reason for making the transformation was to remove the singularity created by the 
term cr^^. The resulting system (pT!) can be analyzed using dynamical systems theory. Namely, 
we can understand the flow near the line of equilibria, thus obtaining insight into the structure 
of solutions of fl2Ul) . We focus on the line of equilibria with x = —1 since it corresponds to 
the attracting x direction. Observe that there exists a center manifold containing the line of 
equilibria. We denote this center manifold by M2d- Note also that points on M2D satisfy the 



constraint a; = — 1 + O(cr^) by center manifold theory 4^. Substituting into fl2T]) . we obtain 

a' = —a^-l - z + 0{a')) 

2 ^ (23) 

z' = a{aji -a + hz + Oia"^) - - ~z + 0{a'^))i). 

Now we can introduce the "inverse" of the time rescaling which lead from fl2Ul) to by 
reducing a factor of a. Thus we are back to the original time coordinate. The resulting system 
is 

1 



a' = --a\-l - ~z + 0{a')) 

z' = crfL- a + bz + 0(0-^) - -o-(-l -z + 0{o-'^))z. 



1 (24) 



2 

System ([23]) has an equilibrium at {a = 0,z = a/b), which has one stable direction and one 
center direction. The branch of the center manifold existing for a > 0, which we denote by 
Mid, is unique. Equation restricted to Mid has the form 

cr' = -V(-l-^ + 0(a)) (25) 

and gives a solution of fl20|) and thus of f|T8|) which is algebraic in t backwards in time. All 
other solutions of ( l20l) grow at least exponentially as t ^ —00. We will denote the special 
solution in the {x, y, z) coordinates by 7^. 

We now look for a special solution of (fT8|) which is a polynomial vector function in t. 
Suppose that 7poiy(i) = ix{t),y{t),z{t)) is a polynomial solution and note that it can only 
exist if x{t) and z{t) are linear in t and y{t) is quadratic in t. We can write 

x{t) = Xq + at, with a some unknown constant. 



16 



Using the first equation in flTSl) . we obtain 

y{t) = (xq + atY — a. 
Substitution into the second equation in (fT8|) results in 

z{t) = (1 - 2a){xo + at). 
Finally, the third equation in f|T8l) gives 

(a — 2a^) = fx + (aa + b{a — 2a'^))t, 
from which we get the conditions 

a — 2c? = Jx 

aa + h{a — 2o?) = 0. 

Solving we get 

a + b _ a{a + b) 

It follows that a polynomial solution exists for Ji = Jiq with 

a{a + b) 
= 2^- 

Note that transformation (|T9|) takes solutions with algebraic behavior as t — > — oo to solutions 
with algebraic behavior as t — > — oo. Consequently, by uniqueness of Mi^), 7^^, = 7poiy We 
denote the value of /x corresponding to flo by /iq. 

Remark 3 We could repeat the construction of the beginning of this section to obtain special 
solutions which are algebraic as t ^ 00. For each p, there exists exactly one such solution. 
Typically the special solutions for t ^ —00 do not match with the special solutions for t —>■ 00. 
The number /iq is an exceptional value of /i for which the two kinds of special solutions are 
equal; the corresponding solution is 7poiy In Section UlIJI we will present numerical evidence 
that 7poiy is the e — limit of a family of canard solutions. 

F. Continuation of SiD.e into the fold region 

In Section IIII CI we proved that all trajectories on Sa^e are attracted to a one dimensional 
manifold Sid^^- Naturally Sio,e and its fiow continuation play a crucial role in the organization 
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of the dynamics of (jS])- In this section, we show that any solution with initial condition on 
SiD,e, when it passes through the region of validity of the {x,y, z) coordinates, must be 0{e) 
close to the special solution 7^. 

As in Section UlIEt we use the coordinate transformation (IT^ . but we apply it to equation 
(IT7|) here. This yields 

x' = (T^^{--a^{x - z)x + (-1 + + -Dix'-^) 
2 a 

a' = -^a\x-z) (26) 

z' = afi + ax + bz — — z)z. 

Subsequently, we apply the time rescaling resulting in the removal of the factor a~^: 

x' = --a^{x - z)x + (-1 + x^ + -Dix^) 
2 a 



'- ^^3/^ ~, (27) 



1 

o = —-o"{x — z 
z' = a{aji + ax + bz — -a{x — z)z 



Notice that setting e = in ([2S]) (resp. ([27])) gives (resp. (I2TD). 
We now note that 



e ^ X ^ z 



^fy ^fy ^fy 

where (x, cr, z) are the coordinates introduced in (fT9l) . Consider a point p G Sxr,^^ close to the 
fold region and let p be the corresponding point in the (x, a, z) coordinates. By (fT6l) . p is close 
to the curve [—p^f?^ ~^^~^)- We fix p > small, which corresponds to the assumption that 
p is close to the fold region. Observe that p is close to (—1, |, — ^^^^)- We investigate the 
solution of fl27j) with initial condition p. 

Consider the following system of equations: 

x' = -\<i^[x - z)x + (-1 + x^ + Dx'jJX?) 

a = --0 (x - z) 

1 (28) 
z = a{ap, + ax + bz — -^'-'^{^ ~ z)z) 

uj = -cr uj(x — z). 
2 ^ ' 

Note that (l28l) has a constant of motion H{a,uj) = ua, and that (127|) is a reduction of fl28l) 
to the invariant hypersurface e = ujct. 

We now argue in a similar way as in section UlI El 
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there exists a center manifold satisfying the constraint a; = — 1 + 0{a'^,uj). We denote 
this manifold by because, even though it is three dimensional in the context of 
fl28l) . it is two dimensional in the variables {x,a,z), i.e. once the reduction uj = ^ has 
been implemented. The manifold can be viewed as a center manifold of the 2D 
surface of equilibria defined by a = and x = —1 + 0{uj). 

restricted to we can rescale the time and find a 2D center manifold, denoted by 
Mfj-), which is a natural generalization of Mi£, to the e > case. That is, center manifold 
theory guarantees the existence of functions ipj{a,uj) such that Mf^ is defined by the 
equations 

X = -1 + 'ipi{a,uj), z = - + 'ijj2{a,uj), (29) 
and Mid is defined by the equations 

x = -l + ^i((T,0), z = ^ + M^,0). (30) 

(Recall that Mid corresponds to the special solution.) The flow on M^^ is given by 

a' = -ia^(-l-^ + 0(^,a)) 

,1 a (31) 

uj' = -auj{-l - + 0{uj, a)) 

and the level curves of the function H = uja are invariant. (This system is the analogue 
of fl25l) in the e = case.) 



We define two sections of the flow of (l28l) : 

S*" = {(x,a,z,p), S-* = {(x,p,z,^), 

where p is the constant used in the definition of the point p. (See Figure [131) We are interested 
in describing the transition by the flow of (|28|) from E*" to We state a few important 

properties of 



1. The manifold is strongly attracting. To see this, recall that is an attracting 
center manifold of the line of equilibria defined byx = — 1, (T = 0. 

2. The manifold Mf^ is also attracting, since it is an attracting center manifold in the 
context of the system obtained by restricting (128|) to and cancelling a factor of a 
(similar construction as leading from fl23|) to (!24l) ). The attraction rate is weaker, but 
since the solutions on Ml^i niove fairly slowly, solutions with initial conditions off Mf^, 
in S*" are very close to Mf^, in E""*. 
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3. Consider a trajectory (f){t) of ( l28l) starting in S™. Since is contained in the invari- 
ant hypersurface ua = e, it follows that it must approach the intersection of ua = e 
with the attracting manifold Mf^j. Since this is a one dimensional curve containing no 
singularities, it must be a trajectory. Hence it follows that approaches a specific 
trajectory on Mf^. 

4. By the definition of and due to the constraint e = auj we have fl = 
{(—1 + p)) P' f + V'2(p, |), |)}, which, together with ( l30l) implies that the distance 
between n Mf^ and S'"'* n Mid is 0(e). Let T be defined by the requirement 
{(t){T)} = n Mf^. It follows that 0(T) converges to 7^ as e ^ 0. 

It follows from points 1. and 2. that the transition map from S*" to defined by the 
flow of (!27|) . mapping a neighborhood of p to a neighborhood of 0(T), is a strong contraction. 
We can conclude this because the restriction of the flow of (!28|) using the constraint e = uoa 
to the flow of ( 1271) amounts to factoring out the neutral direction. 

Remark 4 The manifold as we defined it (center manifold of the surface of equilibria 
given by (7 = and x = — 1 + 0(0;) for (1281) ) is not unique. We add a requirement that 
equals Sa^e in the section S*". This way is the continuation of Sa^e in the (x,(T, 5) 
coordinates. Recall that is strongly attracting, which means that it inherits the stability 
property of Sa.e- The procedure of extending slow manifolds into the flow region is described 
in detail in [3^. 

G. The manifolds M^'^ and M^^^ 

In the preceding section we constructed and analyzed the manifold M^^^ as the continuation 
of Sa,e in the (x, a, z) coordinates. In this section we carry out an analogous procedure, namely 
we find a manifold Mg]^ that extends the unstable slow manifold Sr,e backward in time. Note 
that (!28|) has a equilibria given by the conditions cr = and x = 1 + 0{uj). Let M^Jj be a 
center manifold of this surface of equilibria. This manifold satisfies the condition 

x=l + 0(a^cj). (32) 

The manifold M^'^ is foliated by two dimensional surfaces auj = e, defining invariant sub man- 
ifolds of ([27D. 

We now substitute the constraint (1321) into the (a, z, uj) subsystem of ( l28l) and cancel a 
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factor of cr, analogously to the derivation of fl2^ from fl?Il) . We obtain: 

a' = -^cx\l-z + 0{cx,u)) 

z' = aji + a(l + 0(a^ uo)) + hz- ^aS{l -z + 0{cr, u)) (33) 
u' = -<yuj{l — z + 0{a, u)). 

The system (1331) has an attracting two dimensional invariant manifold M^'^ given by condition 
(1321) and 

z = -- + 0(a,uj). 



The manifold M^'^ defines one dimensional submanifolds of ( l28l) via the constraint e = auo. 
The manifold M^'^, as defined, is not unique. We add the requirement that M^'^ is equal to Sr,e 
in the section S*" (see Remark H]). A very important feature of is that z is contracting 
forwards in time, which means that if the manifold is continued backwards in time, the z 
direction is rapidly expanding. In the e — > limit there exist manifolds and of (l2Til 
that are analogous to M2D and Mir,. 

H. A bound for the MMO behavior 

As mentioned in Section IIIID| MMOs cannot exist if Ji is too large. In other words, the 
region of the the existence of MMOs in the parameter /i must be 0{e) in size. Here we state 
and prove the relevant result. 

Proposition 1 There exists L > and > such that for e < Eq and p, > L no MMOs can 
exist. 

Proof By the analysis of Section llll Fl we know that solutions originating in SiD,e, transformed 
to the variables (x, a, z), satisfy x = —1 + 0{a^, u) and 5 = | + 0{a, u). By assumption ( |T5|) . 
we have z > x for sufficiently small a and w. By possibly decreasing the bound on a, we 
can make sure that z > x is satisfied even when we increase ft (see fl28|) ). We transform 
to the [x, y, z) variables at cr = 5, where 5 is very small and possibly depends on fi. The 
corresponding value oi y is y = 1/5"^ = K. We carry out the remainder of the argument for 
system (fTSjl . leaving the generalization to (fT7|) to the reader. 

We choose a constant Ki satisfying > i^i > so that we can make a "tube" around 
the parabolic cylinder y = x^, with the top and the bottom given hy y = K and y = Ki 
so that the solutions can exit the tube only through the bottom. We define as sides of the 
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tube the parabolic cylinders y = {1 ± a)x'^, where a > is a small constant, and leave the 
tube unbounded in the z direction. For the moment we restrict attention to the part of the 
trajectory for which z < (and x < 0). We now prove that if K and Ki are sufficiently large, 
the trajectories cannot exit the tube through the sides. Note that due to ^ < in flTSl) . we 
obtain x' = ~^ax^ on ?/ = (1 ±a)x^ with \y'\ < \x\ in both cases. If Ki is sufficiently large, the 
sides of the tube are almost vertical, and the (x, y) component of the vector field restricted 
to the sides of the tube is almost horizontal and must be pointing towards the interior of the 
tube. It follows that trajectories can only exit through the bottom, as claimed. By increasing 
ft we can guarantee that z = when the trajectory reaches the bottom of the tube. 

Now note that z will remain non-negative. Also, as long as x remains non-positive, y must 
decrease and z < —fi/b (since z' < fl + bz). However, x cannot reach as long as y is positive. 
Also, since y < Ki we must have x > —y^Ki/ (1 — a). Hence, y' > —{^/Ki/{1 — a) — Ji/h). 
It is now easy to estimate that the time T need for y to reach satisfies 



- v/^i/(l-«)-f 

By taking Ji sufficiently large we can guarantee that by the time y reaches 0, z must become 
at least —hKi/2. By increasing /i we can also make Ki arbitrarily large. Since z will not 
decrease, it follows that x must become at least equal to —bKi/2 before y starts increasing. It 
is well known that for any initial condition (x, y, z) with y < 0, z > 0, and x sufficiently large 
the solution must enter a relaxation loop [34|. Hence, no STO (which requires y to increase 
and X to decrease) can exist. I 

I. Canard solutions and the boundaries of MMO behavior 

In this section we discuss the role of canards in transitions between different MMO regimes. 
We explain in more detail what we mean by MMOs and canards and discuss how transitions 
between different MMO regimes can occur. We identify one type of the possible transitions as 
occurring through a canard and present numerical evidence showing that the final transition 
to spiking must be of this type. Finally we discuss the e — > limit of canards in the (x, y, z) 
variables and show that they correspond to singular canards that are special solutions of f|T8|) . 
An explicit example of such a special solution is the algebraic solution -'^poiy found is Section 
nil E[ We conjecture that the Ji value corresponding to "^poiy is the e = limit of the Ji values 
where the transition between MMOs and spiking occurs. 

We begin by giving a more precise meaning to MMOs and canards. Let V be the in- 
terior of the cube [—5, Sf", where 5 > is some fixed small number. Consider a trajectory 
{x{t), y(t), z(t)) originating near Sa^e outside of V (with the initial value of y > 6). We say that 
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{x{t),y{t), z{t)) has rotation type k if x{t) has k non-degenerate maxima and minima during 
the passage of the trajectory through V. A trajectory originating in Sa^e is a kth secondary 
canard if x{t) has k non-degenerate maxima and minima during the passage through V and 
subsequently continues into Sr^s- A primary canard is a trajectory that passes through V with 
x{t) monotonically increasing (x'(t) > 0) and subsequently continues into S'j.,^. 

Before we discuss MMOs and their transitions, we point out a difference between (1141) and 



the folded node problem |29l . l30l . |31| . Recall the one dimensional invariant manifold SiD,e 
contained in Sa^e- It was shown in Section Till CI that all trajectories on Sa,e are attracted to 
SiD,s- We define an Sm canard as a trajectory originating in Sid^^ and continuing to Sr^^- 
Note that an Sin canard corresponds to an intersection of a ID invariant manifold Sio,e aiid 
2D invariant manifold Sr^e- If such an intersection is transverse with respect to the variation 
of /2 (i.e. in the 4D extended phase space including p,), then SiD,e crosses from one to the 
other side of Sr,e- If "we knew that all the intersections were transverse, we could conclude 
that SiD canards exist only for isolated values of fl. Suppose that 7(t) is a canard solution 
that is not an Sm canard (the initial condition is not contained in SiD,e)- Since 7(t), during 
its passage along S'a,e, becomes exponentially close to SiD,e, it follows that the (minimal) 
separation between SiD,e and Sr,e niust be exponentially small. Consequently, if we knew that 
the distance between Sid^s and Sr^e as a function of p, were not identically equal to and 
varied at least algebraically in fl, we could conclude that canard trajectories can only exist on 
exponentially small (0(e^^/'^), c = const > 0) intervals of the parameter fi. By contrast, in 



the 
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bided node problem, canards exist robustly (for all the values of the control parameter) 
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3li. 



We now consider transitions between different types of MMOs. We will restrict attention 
to MMO periodic orbits of type 1*^, i.e. k small loops followed by a spike. More precisely, we 
consider MMO periodic orbits with rotation type k (see the definition above). If we change ft 
then the rotation type can change in two ways: 

(a) Either a large loop can decrease in size and become fully contained in V, or a loop 

contained in V can grow and be no longer contained in V. 

(b) For some intermediate value of p., x{t) has a degenerate critical point corresponding to 

the collapse of a small loop. 

The following result shows that if a certain non-degeneracy condition holds, then all the 
transitions are of type (a). Moreover, they correspond to transitions through a canard. 

Proposition 2 Suppose there exist two MMO periodic orbits with rotation type k and k — 1, 
respectively, for two different values of p, denoted by pi and p2, respectively. Recall the special 
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solution 7^ = {xfi(t),y^(t),Xp,(t)) and assume that for any ft G [/ii,/i2] ^pit) has only non- 
degenerate critical points. Then, provided that e is sufficiently small, the transition between 
the two MMOs is through a canard. 

Proof Let 71 and 72 denote the MMO with rotation type k and k — 1, respectively. Suppose 
we vary Ji from Jii to /i2. Given a trajectory {x{t),y(t), z(t)) of ([5]), note that a critical point of 
x{t) corresponds to a crossing with the nuUsurface y = x"^ + Dix^. At such a critical point we 
have x" = —y' = —6'^{x — z), so that a degenerate critical point can only occur if x = z. If the 
solution crosses the nullsurface y = x"^ + Dix^ outside of V and close to Sr^e^ then the solution 
must have followed Sr^s for a significant amount of time so that fi + ax + bz = 0{6). Thus 
ax + bz = O^e), by Proposition [H Hence, for e sufficiently small, the critical point must be 
non-degenerate. Now we consider a crossing of the x-nuUsurface inside of V. Any such crossing 
point must also be non-degenerate, due to the assumption on 7^ (by compactness and smooth 
dependence on initial conditions, the intersections of Mfj^ with y = + eDix^ must also be 
non-degenerate). For (x, y) outside of a bounded region we switch to the {x, a, z, u) coordinates 
and observe that, in order to reach the nullsurface, the trajectory would have to follow M2^, 
and therefore be attracted to M[]^, so that z ^ —a/b would have to hold (see Section UlI Gp . 
making x = z impossible since a < —b. We have proved that all the crossing points with the 
y = x^ + Dix^ nullsurface must be non-degenerate. It follows that the only way a transition 
from 7i to 72 can occur is if a point of intersection of the trajectory with the nullsurface 
y = x"^ + Dix^ traces the entire nullsurface until one of the small loops develops to a spike. I 
Although we believe that all the intersections of 7^ with the nullcline y = xP' (for system 
f fTSj) are transverse, we have only partial results in this direction, and we feel that they are 
not suitable for this article and will be presented in future work. Here we state the following 
conjecture. 

Conjecture 1 The least t for which = (for system ( |T8i) ) corresponds to a transverse 
intersection. 

As support for this conjecture we have plotted the value of x — ^ at the first intersection point 
of 7^ with the surface y = x^ m. Figure [HI 

Remark 5 Suppose that conjectureU\is true. Then modifying slightly the proof of Proposition 
IE we could prove that the transition from MMO to spiking is through a primary canard. 

The following result concerns the e — limit of canards in the context of ( JTTl) . We say that 
the special solution 7^ f|T8|) is a singular canard if it continues to M2^. 

Theorem 1 Suppose Ji{e) is a continuous function such that for each {fL,e) satisfying = 
Ji{e) system ( |T7|) has a primary canard. Then, in the e — > limit, system ( |T8l) has a singular 
canard. 
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Proof Recall that in Section fill Fl we showed that the continuation of SiD,e to the fold region 
must be 0{e) close to the special solution 7^. As shown in Section UlI G\ equation fl26l) has 
a 2-d repelling invariant manifold M^^, which is a continuation of Sr^e to the fold region 
(backward in time). In the limit e ^ we obtain a two dimensional invariant manifold of 
( I2U]) . which we denote by Mj^,. A canard solution, as it passes through the fold region, must 
become exponentially close to Mf^, and then connect to M^'^, before going on to Sr,e- In the 
e ^ limit we get a connection from 7^ to Mg^. The /i value corresponding to an e — limit 
of a transition from MMOs to spiking has to be one for which such a connection exists. Note 
that the validity of this argument relies on the convergence of the continuation of SiD,e to the 
special solution 7^ in the variables {x,a,z). The relevant convergence result was established 
in Section UTTFI ■ 

Remark 6 The defining property of a singular canard, namelyy/i continuing to M2D, is non- 
generic and is likely to occur for isolated values of Ji (provided that a suitable transversality 
condition holds). We know one example, with Ji = Jiq = — corresponding to the poly- 

nomial vector solution 7poiy = 7^0 ■ conjecture that fto is the e limit of the locus of 
the transition from MMO to spiking. This conjecture will be presented in more detail at the 
beginning of the next section. 



J. Transition from MMO to spiking — numerical results 



We propose that "jpoiy is the unique singular canard and that it persists as a regular Sid 
canard, that is a connection between M^^ and M^'^- Consequently, we conjecture the following: 



Conjecture 2 There exists a function fi{e) such that for fi = fi{e) system ([5]) has a pri- 
mary SiD canard. The transition from MMOs to spiking occurs within an exponentially small 
(0{e~^^'^), c = const > 0) interval of ft about ft{e). Moreover, 

a{a + b) 

^(0) = /io = ^ — . 

We investigate this conjecture numerically by first considering the transition threshold be- 
tween MMOs and pure spiking in system (ITTl) by computing solution trajectories for decreasing 
values of £ > and /i near fiQ The simulations show that for decreasing values of e > 0, the 
transition between MMOs and pure spiking occurs for decreasing values of fi and suggest that 
/i -> fi{0) ^ 0.1238928363 as e ^ as seen in Figure [H 

In addition, we have found an approximate primary canard numerically. As noted above, 
such a solution must become exponentially close to the attracting manifold Mf^ as it passes 
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through the fold region and then connect to the repeUing manifold M^'^ before going on to Sr,e- 
We exhibit this numerically using the Poincare section option of XPPAUT (Ermentrout02) 
for 6 = 8(1- That is, we choose an initial condition with x = —2 on the singular primary 
canard 7poiy and integrate until reaching the Poincare section x = 0. See Figure dSl This 
solution is expected to become exponentially close to M[^. Next, we approximate M^'^ by 
choosing a range of initial conditions varying in z with x = 5, y = 10, and z near 1.917386 
and integrating backwards in time up to the Poincare section at x = to form a family 
of curves which form an approximation of M2^. As shown in Figure [16], an approximate 
connection between these manifolds is obtained for /2 = .1381943. The approximate primary 
canard solution that corresponds to the intersection of these two manifolds and shows the 
characteristic long term (0(1)) following of a repelling manifold is also displayed in Figure [T71 

Additionally, numerical simulations with varying levels of ft show a transition between 
MMOs and spiking near the value of /i for the approximate primary canard as seen by the 
MMO solution trajectory for = .138 and the spiking solution for fl = .139 in Figure [T71 
Thus, our numerical simulations support the validity of the persistence of the primary canard 
as stated in Conjecture O 

There are many other interesting issues to analyze in this model including the transition 
from small amplitude dynamics to MMOs and the rotation properties of MMOs. For example, 
one could try to obtain an approximation for the number of subthreshold oscillations per spike 
for a given value of /2. We discuss such issues in the Discussion, but a detailed analysis of 
these issues is beyond the scope of this paper. 



IV. DISCUSSION 



In this work, we have analyzed a model introduced by Acker et al. [15| to study synchro- 
nization properties of stellate cells in layer II of the entorhinal cortex (EC). We were interested 
in studying subthreshold oscillations in pyramidal cells of layer V of the EC. Since there was 
no available model of EC layer V pyramidal cells at the time, and since the two cell types 
exhibited very similar electrophysiological characteristics, including the presence of persistent 
sodium and slow potassium currents that are thought to underlie the dynamics of STOs, we 
chose to employ the Acker model to model STOs in layer V pyramids. We note that layer 
V pyramidal cells exhibit many other interesting dynamical properties. They participate in 



high frequency (~ 200 Hz) "ripple" activity associated with hippocampal sharp waves 451] and 



are a site of initiation and propagation of epileptiform activity [46|]. Also, they exhibit graded 
persistent activity, which is associated with memory processes, in which neurons respond to 
consecutive stimuli with graded changes in firing frequency that remain stable after each stim- 
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ulus 43] . In 2006, Fransen et al. introduced a quite detailed multi-compartmental model of 
layer V pyramidal cells that exhibit the graded persistent activity. In future work, we plan 
to adapt our analysis of the simpler model considered here to this more realistic model, and 
analyze the dynamics of the graded persistent activity. 

We began our work with a numerical investigation of the dynamical properties and bifur- 
cation structure of the full model and displayed the role of the persistent sodium and slow 
potassium currents in the generation of subthreshold oscillations. In order to analyze the 6-d 
system, we performed a reduction of dimension analysis to reduce the system to a 3-d system. 
We note that a more thorough reduction of dimensions analysis of a similar model of a layer 
II stellate cells with an h-current and persistent sodium current was performed by Rotstein 
et al. [l2|, and a computational tool for reducing systems possessing multiple scales has also 



been developed 2J]. Rotstein et al. have also analyzed MMOs present in that model 131 ]. 

Just as experimental neuroscientists often modulate or block various currents in order to 
learn about their role by observing the effects of their absence, we investigated the effect of 
modulating various currents in our model. Upon removing or reducing to a very low level 
the persistent sodium current and replacing it with a tonic excitatory input, we witnessed 
very interesting dynamical phenomena including MMOs and seemingly chaotic firing patterns. 
(Future work includes the investigation of whether our analysis and the dynamical mechanism 
involved carry over to the case where a small persistent sodium current is retained.) This 
experiment is unable to be performed presently in the lab since there are no known drugs 
which are able to selectively block persistent sodium current without blocking the standard 
sodium current, and we felt it would be interesting to investigate this regime, especially due to 
the rich dynamical behavior present. The main focus of this work is to explain the dynamical 
mechanisms responsible for the mixed-mode oscillations. 

We have discovered and analyzed a novel dynamical mechanism for the generation of MMOs, 
involving the canard phenomenon as well as the dynamics near a saddle focus equilibrium, 
as in the Shilnikov scenario 48| or the sub critical Hopf scenario 2m, l28|. The techniques we 
employed to analyze the MMOs involve standard techniques from dynamical systems th eory 



such as geometric singular perturbation theory [3J, |49|, |50|] and center manifold theory |44l |. 
For a summary of the analysis, see the last paragraph of section [H 

A main accomplishment of this work is an analytical expression for the parameter at which 
the transition between MMOs and spiking occurs in the singular limit and a characterization of 
the corresponding singular primary canard solution responsible for the transition. Moreover, 
we construct a well developed geometric picture of the dynamics for this new canard-based 
MMO mechanism since we utilize a great deal of the underlying geometric structure available 
including slow manifolds, canards, and the saddle-focus equilibrium with its weak unstable 
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manifold. In addition, the use of the folded-node point from the two time scale model to 
derive our novel three time scale model is new. Also, our analysis, which resulted from the 
investigation of a neuronal model, is applicable to studying MMOs in other systems from 
various contexts that have the same underlying dynamical structure, i.e. systems that can be 
transformed into the form of 

rominent among scenarios for MMOs not involving canards are the Shilnikov scenario 0, 
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5l[ | and the subcritical Hopf scenario 27|,|28||. Both involve a saddle focus equilibrium with a 



two dimensional unstable manifold and a return mechanism of 'soft reinjection', which basically 
means that the return mechanism is given by a typical orientation preserving diffeomorphism. 
Our model system ([5]) has a saddle focus equilibrium with a weak unstable manifold, but 
in addition it has a multiple time scale structure, folded critical manifolds, and a singularity 
similar to a folded node. One of the consequences is that the return mechanism is not soft, but 
has a definite structure determined by the slow manifolds. Moreover, as we have demonstrated 
in the paper, the multiple time scale structure influences subthreshold oscillations by means 
of canard solutions. 

It follows from our analysis that all the orbits leaving the vicinity of the saddle focus 
equilibrium will be attracted to the super-stable slow manifold, which returns to the vicinity 
of the equilibrium. Thus the system can be viewed as close to a homoclinic bifurcation. In this 
paper we have concentrated on the aspects of the problem related to the canard mechanism, 
but we feel that the methods employed to study homoclinic bifurcations, and especially the 
passage near a saddle point, may be applicable and useful for a more detailed investigation of 
our problem. We intend to explore the homoclinic-like features of the problem in future work. 

In future work, we also intend to analyze the rotational properties of the solutions of ([5]). 
For a solution originating in Sa,e or nearby, we define its rotation number to be the number of 
small loops it makes while passing through the fold region, which we define as a sufficiently large 
ball centered at the origin in the coordinates (x, y, z). Note that due to the strong attraction 
of SiD,e 1 all the solutions entering the fold region will have the same rotation number, unless 
the system is very close to a canard. Here we mean a primary canard or a secondary canard, 
as secondary canards correspond to boundary curves, in the (/i, e) parameter space, between 
parameter regions for which solutions have different rotation numbers. We conjecture that for 
a fixed e, there is a sequence of p, values, corresponding to k secondary canards with increasing 
k and to transitions corresponding to an increment by 1 in the rotational number, ending with 
a region corresponding to the maximal number of rotations beyond which only small amplitude 
dynamics exists. Also, we have not analyzed the transition from small amplitude dynamics to 
MMOs. We know from simulations that small amplitude chaos is present, and we conjecture 
that it is related to the loss of normal hyperbolicity of the unstable manifold of the saddle- focus 
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equilibrium, which must occur sufficiently far away from the equilibrium. Also, we believe that 
it is related to the Shilnikov scenario involving incomplete homoclinicity. 

Another interesting question is how to estimate the rotation number for given values of 
and e. Unfortunately, we do not foresee that an analytic formula can be obtained, in contrast 
to related cases in which there are two timescales [30|. A conjecture very strongly supported 
by the numerics is that the rotations happen almost exclusively during the slow drift of the 
solutions along the unstable manifold of of the saddle focus equilibrium. We also believe that 
the flow along this unstable manifold is close to its linear approximation. Once the solutions 
come to a region where x and/or y are sufficiently large, the unstable manifold loses normal 
hyperbolicity and the multiple time scale nature of the system takes over, thus ending the 
rotational phase of the solution. However, in order to produce an approximate formula for 
the number of rotations, one would need to know the distance from the continuation of SiD,e 
to the stable manifold of the saddle focus equihbrium, measured in a specified section such as 
the hyperplane x = 0, as well as the distance between the equilibrium and the region where 
normal hyperbolicity is lost. One would also need the information on the eigenvalues and 
the eigenvectors of the saddle focus equilibrium. No analytic version of the formula can be 
expected since it seems quite unlikely to obtain an analytic formula for either the continuation 
of SiD^ir all the way to the section x = or the stable manifold of the saddle focus equilibrium, 
but a numerical formula seems achievable. 

Note that even in the case of 30|], only a formula asymptotic in e is available, and in fact 
the formula of Szmolyan and Wechselberger is the e = approximation. In our case, we 
can consider an asymptotic formula based on the linear approximation of the flow along the 
unstable manifold, and in the e = approximation this would mean that 7^ would be used 
instead of the continuation of SiD,e- However, again no analytic formula can be expected as no 
analytic expression for 7^ is available, except for the special case where fi = fiQ. A possibility 
would be to obtain an asymptotic expression using center manifold theory, along the ideas of 
Section IIIIEl and we plan to pursue this idea in future work. 
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V. APPENDIX 



A. Model equations 

The full model for the layer V pyramidal cell is given by 

C— = -gNam^h{v - Eno) - gKn^{v - Ek) - gL{v - El) 

-9NapP{v - ENa) " QKsWiy - Ek) + lapp 

dm moo{v) — m 

dt Trn{v) 

dn ^ n^{v) - n ^^^^ 

dt Tn{v) 

dh hoo{v) — h 

dt Th{v) 

dp _ Poo{v) -P 
dt 'Tp{v) 
dw Woq{v) — w 

The default parameters for the system are: C = 1.5, gNap = -21, gNa = 52, gx = 11, 
gi = .1, gxs = 2.0, = 90, E^a = 55, Ek = —90, and El = —54. The graphs of the 
voltage dependent activation and inactivation functions as well as the voltage dependent time 
constants are given in Figures [1] and [21 and their equations are the following: 

m^{v) = arr,{v) / {am{v) + Pm{v)) where 

a^{v) = -.l{v + 23)/(exp(-.l(f + 23)) - 1) and (3m{v) = 4exp(-(f + 48)/18). 
noo(f) = an{v)/{an{v) + Pn{v)) where 

a„(t.) = -.01{v + 27)/(exp(-.l(t; + 27)) - 1) and Pn{v) = .125exp(-(t; + 37)/80). 
hoo{v) = ah{v)/{ah{v) + (3h{v)) where 

a^,{v) = .07exp{-{v + 37)/20) and (3h{v) = l/(exp(-.l(t; + 7)) + 1). 

= 1/(1 + exp(-(t; + 38)/6.5)). 
w,^{v) = 1/(1 + exp(-(t; + 35)/6.5)). 

The equations for the voltage dependent time constant for n is: 

Tr,{v) = l/{an{v) + f3n{v)). 
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VI. FIGURES 
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FIG. 1: Steady state activation and inactivation curves for the gating variables. 
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FIG. 2: Voltage dependent time scale curves for the gating variables. Note: Tw{v) = 90. 
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FIG. 3: Bifurcation diagram for the full system ([I])-© as function of /( 



33 



w 

0.1 




0.02 I 1 1 1 1 1 1 1 1 

-80 -60 -40 -20 20 40 60 

V 

FIG. 4: Solution trajectories of the full system ([H)-(l2]) with lapp = -9 projected onto the w-v plane 
showing bi-stability of spiking solutions and damped subthreshold oscillations. 
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FIG. 5: Currents in the full system ([I])-® with lapp = .9 during damped subthreshold oscillations 
showing the interaction of I Nap and Iks- 

V 

-40 



-45 
-50 
-55 
-60 



1 2 3 4 5 

lapp 

FIG. 6: Bifurcation diagram for the full system (HI)-© with I^a removed as function of lapp- Note 
the change in the location and criticality of the Hopf bifurcation point versus Figure [3l 
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FIG. 9: Bifurcation diagram for the 3-d system ([3]) with the persistent sodium current removed. A 
sequence of period doubhng bifurcations occur starting around lapp = 16.97. 
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FIG. 10: MMO trajectory for Iapp=17.1 in ([3]). Trajectory in blue. Fold curve and folded singularity 
in green. Fixed point and its 1-d stable eigenvector in red. Folded surface is the v-nullsurface, and 
the shaded plane is the unstable eigenspace of the fixed point. (Note that the trajectory pierces the 
v-nullsurface near its fold. This is typical for passage near a folded node or folded saddle node type 
of singularity.) 
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FIG. 11: Caricature of the folded slow manifold. 



39 




-0 . 3 



FIG. 12: SiD plotted along with a spiking solution of ([5]) with £ = .1 that tracks Sid quite closely 
as well as an MMO solution for the larger value of e = £d = .2764436178. /i = .02839. 
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FIG. 15: Values of fl corresponding to the transition between MMOs and pure spiking in system ([5]) 
for decreasing values of epsilon. 
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FIG. 16: Approximation of the primary canard formed by the connection of approximations of Mfj^ 
and M^'^ for fl = .1381943 in system ^ with the default value of e = .2764436178. The family of 
curves approximating the manifold is obtained by taking a sequence of initial conditions with 
X = 5, y = 10 and z near 1.917386 and evolving backwards in time until reaching the x = plane. 
The single curve approximating Mf^ is generated by evolving in forwards time until reaching the 
X = plane with initial conditions x = —2, y = 3.726472, z = —.9058866 on the special solution 
7poiy(i) with y and z as prescribed in Section [III E[ 
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FIG. 17: Top: Approximate primary canard solution of system (|17p with p, = 0.1381943 near the 
transition between mixed-mode oscihations and spiking. Middle: Approximate primary canard so- 
lution of system (jlTD with fl = 0.1381943 near the transition between mixed-mode oscillations and 
spiking. Bottom: Periodically spiking solutions "l^ihibited by system (jl7p with jl = 0.139. 
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